Joint Density-Functional Theory of the Electrode-Electrolyte Interface: Application 
to Fixed Electrode Potentials, Interfacial Capacitances, and Potentials of Zero Charge 
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This work explores the use of joint density-functional theory, a new form of density-functional 
theory for the ab initio description of electronic systems in thermodynamic equilibrium with a liq- 
uid environment, to describe electrochemical systems. After reviewing the physics of the underlying 
fundamental electrochemical concepts, we identify the mapping between commonly measured elec- 
trochemical observables and microscopically computable quantities within an, in principle, exact 
theoretical framework. We then introduce a simple, computationally efficient approximate func- 
tional which we find to be quite successful in capturing a priori basic electrochemical phenomena, 
including the capacitive Stern and diffusive Gouy-Chapman regions in the electrochemical double 
layer, quantitative values for interfacial capacitance, and electrochemical potentials of zero charge 
for a series of metals. We explore surface charging with applied potential and are able to place 
our ab initio results directly on the scale associated with the Standard Hydrogen Electrode (SHE). 
Finally, we provide explicit details for implementation within standard density-functional theory 
software packages at negligible computational cost over standard calculations carried out within 
vacuum environments. 

PACS numbers: 71.15.Mb,82.45.-h,82.45.Gj,73. 
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I. INTRODUCTION 



Ab initio calculations have shed light on many ques- 
tions in physics, chemistry, and materials science, includ- 
ing chemical reactions in solution^'^ and at surfaces^^r— 
However, first principles calculations have offered less in- 
sight in the complex and multi-faceted field of electro- 
chemistry, despite the potential scientific and technolog- 
ical impact of advances in this field. Because the fun- 
damental microscopic mechanisms involved in oxidation 
and reduction at electrode surfaces are often unknown 
and are difficult to determine experimentally,'' rich sci- 
entific opportunities are available for theoretical study. 
From a technological perspective, practicable first prin- 
ciples calculations could become a vital tool to direct the 
experimental search for better catalysts with significant 
potential societal impact: as just one example, economi- 
cally viable replacement of gasoline powered engines with 
fuel cells in personal transport systems requires systems 
operating at a cost of $35/kW, whereas the current cost is 
S294/kW,''. due mostly to the expense of platinum-based 
catalyst materials. 

The primary challenge which distinguishes theoretical 
study of electrochemical systems is that including the liq- 
uid electrolyte, which critically influences the functioning 
of the electrochemical cell, requires detailed thermody- 
namic sampling of all possible internal molecular con- 
figurations of the fiuid. Such critical influences include 
(a) screening of charged systems, (b) establishment of an 
absolute potential reference for oxidation and reduction 
potentials, and (c) voltage-dependence of fundamental 
microscopic processes, including the nature of reaction 
pathways and transition states. While there have been 
attempts at the full ab initio molecular dynamics ap- 
proach to this challenge fi'^'S such calculations are nec- 



essarily of the heroic type, require tremendous computa- 
tional resources, and do not lend themselves to system- 
atic studies of multiple reactions within a series of many 
candidate systems. Such studies require development of 
an alternate approach to first-principles study of electro- 
chemistry. 



A. Previous approaches 

One response to the aforementioned challenges is to 
avoid the issue and lessen the computational cost either 
by forgoing electronic structure calculation entirely or by 
neglecting the thermodynamic sampling of the environ- 
ment. Some studies have employed classical molecular 
dynamics with interatomic potentialsj^.U however, such 
semi-empirical techniques often perform poorly when de- 
scribing chemical reactions involving electron-transfer, 
which are central to oxidation and reduction reactions. 

The latter approach - single configuration ab initio cal- 
culations - neglects key phenomena associated with the 
presence of an electrolyte liquid in equilibrium. The most 
direct single configuration ab initio approach pursued to 
date is to study the relevant reactions on a surface in 
vacuum and to study trends and correlations with the 
behavior in electrochemical systemsJ^ii^ Some of these 
studies are done in a constant charge or constant poten- 
tial ensemblei^ to allow variation of the applied electrode 
potential. This approach, however, does not include criti- 
cal physical effects of the electrolyte such as the dielectric 
response of the liquid environment and the presence of 
high concentrations of ions in the supporting electrolyte. 
In response, an intermediate approach is to include a 
layer or few layers of explicit water molecules into the 
calculationJ^^i^ 
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Such an approach is problematic for a number of rea- 
sons. First, actual electrochemical systems can have 
rather long ionic screening lengths (30 A for an ionic con- 
centration of 0.01 M), which would require large amounts 
of explicit water. Second, simulation of the actual effects 
of dipolar and ionic screening in the fluid requires ex- 
tensive sampling of phase space, corresponding to very 
long run times. Indeed, in some references, only one 
layer of frozen water without thermal or time sampling 
is includedJ^ Moreover, as most reactions of interest oc- 
cur at potentials away from the potential of zero charge, 
such calculations must include a net charge, which can 
be problematic in typical solid-state periodic supercell 
calculations. One may compensate for this charge with 
a uniform charged background extending throughout the 
unit cell, both the liquid and the solid regions, but 
this distribution does not reflect the electrochemical re- 
ality. Other methods include an explicit reference elec- 
trode with a corresponding negative surface charge to 
keep the unit cell neutral,™ but this requires a some- 
what arbitrary choice of where to place the compensating 
electrode and may not lead to realistic potential profiles. 
More recently, modeling the electrolyte by a layer of ex- 
plicit hydrogen atoms was shown to provide a source of 
electrons for charged surface calculations while keeping 
the unit cell neutral. Again, however, this approach 
requires either judicious choice of the locations of the 
corresponding protons which make up the corresponding 
reference electrode or computationally intensive thermo- 
dynamic sampling. 



Another broad approach constructs an approximate 
a posteriori continuum model^i for both the dielectric 
response of the water molecules and the Debye screen- 
ing effects of the ions and performs ab initio calculations 
where the electrostatic potential is determined by solv- 
ing Poisson-Boltzmann-like equations. ^^"^"'^ Explicit in- 
clusion of a few layers of explicit water molecules and 
ionic species within the ab initio calculations can fur- 
ther enhance the reliability of this approach without dra- 
matic additional computational cost. While including ex- 
plicitly the most recognized physical effects of the elec- 
trolyte, such Poisson-Boltzmann-like approaches do not 
arise from an exact underlying theory. Thus, they may 
disregard physically relevant effects, such as the non- 
locality and non-linearity of the dielectric response of liq- 
uid water and the surface tension associated with forma- 
tion of the liquid-solid interface. We note, for instance, 
that a typical electrochemical field strength would be a 
0.1 V drop over a double layer width of 3 Angstroms, or 
300 MV/m, a field at which the bulk dielectric constant of 
water is reduced by about one-third, strongly indicating 
that non-linear dielectric saturation effects are present 
in actual electrochemical systems, particularly near the 
liquid-solid interface, and ultimately should be captured 
naturally for an ab initio theory to be truly predictive 
and reliable. 



B. Joint density-functional theory approach 

This work begins by placing the aforementioned modi- 
fied Poisson-Boltzmann approaches on a firm theoretical 
footing within an, in principle, exact density-functional 
theory formalism, and then describes the path to in- 
cluding all of the aforementioned effects in a fully rig- 
orous ab initio density functional. The work then goes 
on to elucidate the fundamental physics underlying elec- 
trochemistry and provide techniques for computation 
of fundamental electrochemical quantities from a for- 
mal perspective. The work then shifts focus and intro- 
duces an extremely simplified functional for initial explo- 
ration of the potential of our overall approach for prac- 
tical calculations. The equations which result at this 
high level of simplification resemble those introduced by 
others^'^'?^ from an a posteriori perspective, thus putting 
those works on a firmer theoretical footing and showing 
them in context as approximate versions of a rigorous 
underlying approach. We then work within this simpli- 
fied framework to explore - in more depth than previ- 
ously in the literature - fundamental physical effects in 
electrochemistry, including the microscopic behavior of 
the electrostatic potential near an electrode surface, the 
structure of the electrochemical double layer, differen- 
tial capacitances, and potentials of zero charge across a 
series of metals. The encouraging results which we ob- 
tain even with this highly simplified functional indicate 
that the overall framework is sound for the exploration 
of physical electrochemical phenomena and strongly sug- 
gests that the more accurate functional under present 
development's, will yield accurate, fully ab initio results. 

Section II begins by laying out our theoretical frame- 
work. Section III describes connections between exper- 
imental electrochemical observables and microscopic ab 
initio computables. Section IV introduces a simple ap- 
proximate functional which offers a computationally effi- 
cient means of bridging connections to experimental elec- 
trochemistry. Section V provides specific details about 
electronic structure calculations of transition metal sur- 
faces. Finally, Section VI presents electrochemical results 
for those metallic surfaces obtained with our simplified 
functional and Section VII concludes the paper. The 
appendices include technical information regarding im- 
plementation of our functional within a pseudopotential 
framework. 



II. THEORETICAL FRAMEWORK 

As described in the Introduction, much of the challenge 
in performing realistic ab initio electrochemistry calcula- 
tions comes not only from the need to include explicitly 
the atoms composing the environment but also from the 
need to perform thermodynamic averaging over the loca- 
tions of those atoms. Recently, however, it was proved 
rigorously that one can compute exact free-energies by 
including the environment in a joint density-functional 
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theory frameworkf2Si21 Specifically, this previous work 
shows that the free energy A of an explicit quantum me- 
chanical system with its nuclei at fixed locations while in 
thermodynamic equilibrium with a liquid environment 
(including full quantum mechanical treatment of the en- 
vironment electrons and nuclei), can be obtained by the 
following variational principle;^ 



A= min {G[n{r),{N4r)},V{r)] 

nir).,{N^{r)} 



(frV{r)n{r)} 



(1) 



where G[n{r), {Na{r)}, V{r)] is a universal functional of 
the electron density of the explicit system n(r), the den- 
sities of the nuclei of the various atomic species in the en- 
vironment {Na{r)}, and the electrostatic potential from 
the nuclei of the explicit system V{r). The functional 
G[n{r),{Na{r)},V{r)] is universal in the sense that it 
depends only on the nature of the environment and that 
its dependence on the explicit system is only through the 
electrostatic potential of the nuclei included in V{r) and 
the electron density of the explicit system n{r). With this 
functional dependence established, one can then separate 
the functional into large, known portions and a smaller 
coupling term ultimately to be approximated;^ 



G[nir),{N^{r)},Vir)] 



AA[n{r),{N„{r)},V{r)] (2) 



where Axsinir)] and i}iq[{Na{r)}] are, respectively, the 
standard universal Kohn-Sham electron-density func- 
tional of the explicit solute system in isolation (in- 
cluding its nuclei and their interaction with its elec- 
trons) and the "classical" density- functional for the liq- 
uid solvent environment in isolation. The remainder, 
AA[n{r),{Na{r)},V{r)] is then the coupling term be- 
tween the solute and solvent. 

For Aif5[n(r)], one can employ any of the popular ap- 
proximations to electronic density functional theory such 
as the local-density approximation (LDA), or more so- 
phisticated functionals such as the generalized-gradient 
approximation (GGA).^^ On the other hand, functionals 
fliq[{Na{r)}] for liquid solvents such as water are gener- 
ally less-well developed, though the field has progressed 
significantly over the past few years. For example, one 
recent, numerically efficient functional for liquid water 
reproduces many of the important factors determining 
the interaction between the liquid and a solute, including 
the linear and nonZmear non-local dielectric response, the 
experimental site-site correlation functions, the surface 
tension, the bulk modulus of the liquid and the variation 
of this modulus with pressure, the density of the liquid 
and the vapor phase, and liquid-vapor coexistence^. A 
framework employing such a functional would be more re- 
liable than the modified Poisson-Boltzmann approaches 
available to date, which do not incorporate any of these 
effects except for the linear local dielectric response ap- 
propriate to macroscopic fields. Inclusion of the densi- 
ties of any ions in the electrolyte environment among 



the {Na{r)} is a natural way to include their effects into 
^iq[{Na{r)}] and provide ionic screening into the overall 
framework. 

Finally, developing approximate forms for the cou- 
pling AA[n{r),{Na{r)},V{r)] in ([2]) remains an open 
area of research. In an early attempt, Petrosyan and 
co-workers?- employed a simplified iliq[{Na{r)}] using a 
single density field N{r) to describe the fluid. In that 
preliminary work, because such an N{r) gives no explicit 
sense of the orientation of the liquid molecules, the ten- 
dency of these molecules to orient and screen long-range 
electric fields was included a posteriori into a simplified 
linear (but nonlocal) response function. In a more com- 
plete framework with explicit distributions for the oxygen 
and hydrogen sites among the {Na{r)} the full non-local 
and non-linear dielectric response can be handled com- 
pletely a priori^ Beyond long-range screening effects, 
the coupling AA[n{r), {Na{r)},V{r)] must also include 
effects from direct contact between the solvent molecules 
and the solute electrons. Because the overlap between the 
molecular and electron densities is small, the lowest-order 
coupling, very similar to the "molecular" pseudopoten- 
tials of the type introduced by Kim et al.y^ would be a 
reasonable starting point. Using such a pseudopotential 
approach (with only the densities of the oxygen atoms 
of the water molecules), Petrosyan and coworkers?^ ob- 
tained good agreement (2 kcal/mole) with experimental 
solvation energies, without any fitting of parameters to 
solvation data. Combining a coupling functional AA sim- 
ilar to that of Petrosyan and coworkers with more explicit 
functionals V,iq[{Na{r)}] for the liquid^^ and standard 
electron density functionals AKsin-ir)] for the electrons, 
is thus a quite promising pathway to highly accurate ab 
initio description of systems in equilibrium with an elec- 
trolyte environment. 



III. CONNECTIONS TO ELECTROCHEMISTRY 

Turning now to the topic of electrochemistry, we 
present a general theoretical framework to relate the re- 
sults of ab initio calculations to experimentally measur- 
able quantities, beginning with a brief review of the elec- 
trochemical concepts. 



A. Electrochemical potential 

In the electrochemical literature, the electrochemical 
potential fl of the electrons in a given electrode is de- 
fined as the energy required to move electrons from a 
reference reservoir to the working electrode. This po- 
tential is often conceptualized as a sum of two terms, 
A = lJ-int~F^, where fimt is the purely chemical potential 
(due to concentration gradient and temperature, chemical 
bonding, etc.), $ is the external, macroscopic electro- 
static potential, and F is Faraday's constant. (Note that 
F = Nac has the numerical value of unity in atomic 
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units, where chemical potentials are measured per par- 
ticle rather than per mole.) In the physics literature, 
this definition for Jl (when measured per particle) corre- 
sponds precisely to the "chemical potential for electrons," 
which appears for instance in the Fermi occupancy func- 
tion [e("-^)/'=«^ + 1]"^ 



B. Electrode potential 

In a simple, two-electrode electrochemical cell, the 
driving force for chemical reactions occurring at the elec- 
trode surface is a voltage applied between the reference 
electrode and working electrode. In the electrochemical 
literature, this voltage is known as the electrode potential 
£, defined as the electromotive force applied to a cell con- 
sisting of a working electrode and a reference electrode. 
In atomic units (where the charge of an electron is unity) , 
the electrode potential is thus equivalent to the energy 
(per fundamental charge e) supplied to transfer charge 
(generally in the form of electrons) from the reference 
to the working electrode, assuming no dissipative losses. 
Under conditions where diffusion of molecules and reac- 
tions occurring in the solution are minimal, this energy 
is completely transferred to the electrons in the system, 
causing a corresponding change in the electrochemical 
potential of the electrons in the working electrode. 

An idealized two-terminal electrochemical cell con- 
trols the chemical potential of a working electrode /i^^^ 
through application of an electrode potential £ (voltage) 
between it and a reference electrode of known chemical 
potential p,^^^ (See Figure [TJa)). With the application 
of the electrode potential £, the energy cost to the elec- 
trochemical cell, under reversible (lossless) conditions, to 
move a single electron from the reference electrode to the 
working electrode is dll = — /l'-^'' -I- jl^^^ + £. Here, the 
electrode potential appears with a positive sign, because 
to move a negative charge from the negative to positive 
terminal requires a net investment of energy, and thus 
cost to the electrochemical cell, against the source of the 
potential £. Under equilibrium conditions, we must have 
dU = 0, so that £ = /i(^) - p,^-^K 

As Section IV shows, the electrostatic model which 
we employ for ionic screening in this work establishes 
a fixed reference such that the microscopic electron po- 
tential cf) (the Coulomb potential energy of an electron 
at a given point) is zero deep in the liquid environ- 
ment far from the electrode (See Figure [Ijb)) - imply- 
ing that the macroscopic electrostatic potential $ (which 
differs in overall sign from <j>) there is also zero. A con- 
venient working electrode thus corresponds to electrons 
solvated deep in the fluid, which will have electrochem- 
ical potential //^-^^ = — = where /z^^^ cor- 
responds to the solvation energy of an electron in the 
liquid. Referring the scale of the electrochemical poten- 
tial to such solvated electrons (so that = 0), we 
then have p,^^^ = 0, so that £ = —p,^^K In sum, the 
opposite of the electronic chemical potential in our ab 
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FIG. 1: (a) Schematic of an electrochemical cell. The 

working electrode is explicitly modeled while the 
reference electrode is fixed at zero, (b) Relationship 

between the microscopic electron potential {4>{z)) 
(averaged over the directions parallel to the surface), 
electrochemical potential, and applied potential for a Pt 
(111) surface. The large variations in potential to the 
left of zpt correspond to the electrons and ionic cores 
comprising the metal while the decay into the fiuid 
region is visible to the right of zpt. 



initio calculations corresponds precisely to the electrode 
potential relative to solvated electrons. In practice, the 
choice of approximate density functionals ^liq[{N a{r)}] 
and /S.A[n{r),{Na{r)},V{r)\ sets the value of the elec- 
tron solvation energy; each model fluid corresponds to a 
different reference electrode of solvated electrons. Sec- 
tion VII demonstrates the establishment of the electro- 
chemical potential of such a model reference electrode 
relative to the standard hydrogen electrode (SHE). 



C. Potential of zero charge (PZC) and differential 
capacitance. 

For any given working electrode, a specific number of 
electrons, and thus electronic chemical potential /i, is re- 
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quired to keep the system electrically neutral. The cor- 
responding electrode potential {£ — —fl) is known as 
the potential of zero charge. Adsorbed ions from the 
electrolyte or other contaminants on the electrode sur- 
face create uncertainty in the experimental determina- 
tion of the potential of zero charge. One advantage to ab 
initio calculation is the ability to separate the contribu- 
tion due to adsorbed species from the contribution of the 
electrochemical double layer, the latter being defined as 
the potential of zero free charge (PZFC). Experimentally, 
only the potential of zero total charge (PZTC), which in- 
cludes the effects of surface coverage, may be measured 
directly, and the potential of zero free charge can only 
be inferred;^ Ab initio approaches such as ours allow for 
the possibility of controlled addition of adsorbed species 
and direct study of these issues. 

At other values of the electrode potential £, the system 
develops a charge per unit surface area a = Q/A. From 
the relationship between these two quantities <j{£), one 
can then determine the differential capacitance per unit 
area C = ^ . The total differential capacitance of a metal 
is determined by both the density of states of the metal 
surface CpQg , also known as the quantum capacitance)^ 
and the capacitance associated with the fluid Cjj. These 
capacitances act in series, so that full differential capac- 
itance is given by 

=Cgl +Cj^[^g. (3) 

In typical systems, C'pQg ^ 100 — lOQQjiF / cm? is 
larger than the fluid capacitance (typically ^ 15 — 
100/iF/cm^), so when the two are placed in series, the 
fluid capacitance dominates. 

The fluid capacitance may be further decomposed into 
two capacitors acting in series, 

qi=c^i + c-i, (4) 

as in the Gouy-Chapman-Stern model for the electro- 
chemical double layer;^"— The surface charge on the 
electrode and the first layer of oppositely charged ions 
behave like a parallel plate capacitor with distance A be- 
tween the plates. A indicates the distance from the elec- 
trode surface to the first layer of ions - called the outer 
Helmholtz layer for non-adsorbing electrolytes. The ca- 
pacitance per unit area for this simple model is Ca = ^ , 
analogous to the Helmholtz capacitance. For a gap size 
A ^ 0.5 A, this model leads to a "gap" capacitance of 
about 20 ^F/cw?. Additional capacitance arises from 
the diffuse ions in the liquid, where the model for this ca- 
pacitance = eheoKCOsh(ll^) is also well-known from 
the electrochemistry literaturci^ In the limit where most 
of the voltage drop is found in the outer Helmholtz layer 
((/)(A) ~ fcsT), this expression reduces to a constant 
value which depends only on the concentration of ions 
in the electrolyte and the bulk dielectric constant of the 
fluid Ck = For water with a 1.0 M ionic con- 

centration, the "ion" capacitance is C„ =240 fiF/crr?, an 



order of magnitude larger than the "gap" capacitance. At 
this high ionic concentration, the "gap" (Helmholtz) ca- 
pacitance dominates not only the fluid capacitance, but 
also the total capacitance. For lower concentrations of 
ions, the magnitude of the "ion" capacitance becomes 
more comparable to the "gap" capacitance and voltage- 
dependent nonlinear effects in the fluid could become im- 
portant. 



D. Cyclic voltammetry 

A powerful technique for electrochemical analysis is 
the cyclic voltammogram, in which current is measured 
as a function of voltage swept cyclically at a constant 
rate. Such data yield detailed information about electron 
transfer in complicated electrode reactions, with sharp 
peaks corresponding to oxidation or reduction potentials 
for chemical reactions taking place at the electrode sur- 
face. Because current is a time-varying quantity and 
density-functional theory does not include information 
about time dependence and reaction rates, careful rea- 
soning must be employed to compare ab initio calcula- 
tions to experimental current-potential curves. Previous 
work has correlated surface coverage of adsorbed hydro- 
gen with current in order to predict cyclic voltammo- 
grams for hydrogen evolution on platinum electrodes. — 
This simple model for a cyclic voltammogram is intrinsi- 
cally limited at a full monolayer of hydrogen adsorption, 
rather than by the more realistic presence of mass trans- 
port and diffusion effects, but nonetheless provides useful 
comparisons to experimental data. 

Using a similar approach, our framework gives the pre- 
dicted current density J directly through the chain rule 
as 

at at at 

where K = is the voltage sweep rate, and C{E) is the 
differential capacitance per unit area at electrode poten- 
tial £, as defined above. For the bare metal surfaces with 
no adsorbates studied in Section VII of this work, only 
the double layer region structure is visible, but the tech- 
nique may be generalized to study chemical reactions at 
the electrode surface. The current density curve is sim- 
ply proportional to the differential capacitance per unit 
area C as long as the state of the system varies adiabat- 
ically and the voltage sweep rate is significantly slower 
than the reaction rate. In the adiabatic limit, features in 
the charge-potential curves calculated for reaction inter- 
mediates and transition states can be compared directly 
with peaks in cyclic voltammograms to predict oxidation 
and reduction potentials from first principles. 
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IV. IMPLICIT SOLVENT MODELS 

For computational expediency and to explore the per- 
formance of the overall framework for quantities of elec- 
trochemical interest, we now introduce a highly approx- 
imate functional. Despite its simplicity, we find that the 
model below leads to very promising results for a number 
of physical quantities of direct interest in electrochemical 
systems. The first step in this approximation is to min- 
imize with respect to the liquid nuclear density fields in 
the fully rigorous functionaP^ so that Eq. ([T|) becomes 

A = miiv{AKsHr),{Zi,Ri}] 

n{r) 



AA[n{r),{Zi,Ri}]), 



(5) 



with the effects of the liquid environment all appearing 
in the new term 

l\A[n{r),{Zi,Ri}] = min (fliJ^^W] 

Nc,{r) 

+ AA[n(r),iV„(r),{Z,,i?z}]), (6) 

where Zi and Ri are the charges and positions of the sur- 
face nuclei (and those of any explicitly included adsorbed 
species). This minimization process leaves a functional 
in terms of only the properties of the explicit system and 
incorporates all of the solvent effects implicitly. Up to 
this point, this theory is in principle exact, although the 
exact form of AA[n(r), {Z/, i?/}] is unknown. For prac- 
tical calculations this functional must be approximated 
in a way which captures the underlying physics with suf- 
ficient accuracy. 



A. Approximate functional 

In this initial work, we assume that the important in- 
teractions between the solvent environment and the ex- 
plicit solute electronic system are all electrostatic in na- 
ture. Our rationale for this choice is the fact that most 
electrochemical processes are driven by (a) the surface 
charge on the electrode and the screening due to the di- 
electric response of the liquid solvent and (b) the rear- 
rangement of ions in the supporting electrolyte. To in- 
corporate these effects, we calculate the electron poten- 
tial 4){r) (the Coulomb potential energy of an electron 
at a given point, which equals — e times the electrostatic 
potential) due to the electronic and atomic core charges 
of the electrode and couple this potential to a spatially 
local and linear description of the liquid electrolyte envi- 
ronment, yielding 

A[n{r),(l){r)] ^ ATxc[n{r)] 

+ j <fr{^{r) {n{r)-N{r,{Zi,Ri})) 

g'J^-(0M)n, (7) 



where ATxcl^ir)] is the Kohn-Sham single-particle ki- 
netic plus exchange correlation energy, n(r) is the full 
electron density of the explicit system (including both 
core and valence electrons), and iV(r, {i?/, Zj}) is the nu- 
clear particle density of the explicit solute system with 
nuclei of atomic number Zj at positions i?/, eb is the 
bulk dielectric constant of the solvent, and e(r) and K(r) 
are local values, respectively, of the dielectric constant 
and the inverse Debye screening length due to the pres- 
ence of ions in the fluid. We emphasize that, despite 
the compact notation in ([7]), in practice we employ stan- 
dard Kohn-Sham orbitals to capture the kinetic energy 
and, as the appendices detail, we employ atomic pseu- 
dopotentials rather than direct nuclear potentials, so that 
A^(r, {i?/, Z/}) does not consist in practice of a set of 
Dirac (5-functions. 

To determine local values of the quantities e(r) and 
K(r) above, we relate them directly to the local average 
density of the solvent Niq{r) as 



e(r) 



(8) 



where Nb and are, respectively, the bulk liquid number 
density (molecules per unit volume) and the bulk dielec- 
tric constant, and = ^^j. ^ ■ NiZ^ is the square 
of the inverse Debye screening length in the bulk fiuid, 
where Zi and Ni = CiNa are the valences and number 
densities of the various ionic species. Finally, our model 
for the local liquid density depends on the full solute 
electron density n(r) at each point through the relation 



Niqin) = ^erfc 



In {n/na) 
%/27 



(9) 



a form which varies smoothly (with transition width 7) 
from the bulk liquid density Nb in the bulk solvent where 
the electron density from the explicit system is less than 
a transition value no to zero inside the cavity region as- 
sociated with the solute, defined as those points where 
n{r) > hq. This form for Niq(n) reproduces solvation en- 
ergies of small molecules in water without ionic screening 
to within 2 kcal/mol,^^ when the parameters in Eq. (|9]) 
have values 7 = 0.6 and no = 4.73 x 10^"^ A^'^. 

The stationary point of the functional in Eq. ([7]) de- 
termines the physical state of the system and is actu- 
ally a saddle point which is a minimum with respect 
to changes in n(r) (or, equivalently, the Kohn-Sham or- 
bitals) and a maximum with respect to changes in (f>{r). 
Setting to zero the variation of Eq. ([7]) with respect to the 
single-particle orbitals generates the usual Kohn-Sham, 
Schrodinger-like, single-particle equations with (/>(r) re- 
placing the Hartree and nuclear potentials and results in 
the modified Poisson-Boltzmann equation. 



V-(e(r)V0(r))-e6K2(r)0(r) 
-.-An {n{r)- N{r,{Ri,Zi})). 



(10) 
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FIG. 2: Microscopic and model quantities for Pt(lll) 
surface in equilibrium with electrolyte: (a) Pt atoms 
(white), electron density n(r) (green), and fluid density 
Niq{r) (blue) in a slice passing from surface (left) into 
the fluid (right) with zpt — 5.95 A indicating the end 
of the metal, (b) dielectric constant (e(z)) and screening 
length (k~^(z)) (averaged over the planes parallel to the 

surface) for ionic concentrations of 1.0 M and 0.1 M 
along a line passing from surface into the fluid. Position 
z — zpt measures distance from the end of the metal 
slab. (See Sections V and VI.) 



Self-consistent solution of this modified Poisson- 
Boltzmann equation for (j){r) along with solution of the 
corresponding traditional Kohn-Sham equations defines 
the final equilibrium state of the system. 

Figure [2] illustrates the various concepts in this model 
using actual results from a calculation of the Pt(lll) sur- 
face, described in Sections V and VI. Figure [2Ia) shows 
the electron n{r) and liquid Niq(r) densities in a slice 
through the system which passes through the metal (left, 
z < zpt ) and the fluid (right, z > zpt). We define the 
end of the metal surface ^^j^-^gi^^l covalent radius 



of the last row of metal atoms (zpt = 5.95 A). The ionic 
cores and the itinerant valence electrons in the metal are 
visible, as well as the gap between the surface and the 
bulk of the fluid. As shown in Figure [Ifb), the local 
functions for the dielectric constant e(r) and the inverse 
Debye screening length K{r) respect the correct physical 
limiting values: eb and Kb in the bulk solvent and e = 1 
and K = within the surface. The rapid increase in di- 
electric constant for A< z — zpt < 1 A corresponds 
to the appearance of fluid on the right side and results 
in the localization of significant charge from the fiuid at 
this location. The inverse screening length k depends on 
the concentration of ions in the electrolyte through the 
bulk liquid value Kb- Figure [IJb) shows screening length 
as a function of distance from the metal surface for both 
0.1 and 1.0 molar bulk ionic concentrations. The large 
screening lengths at positions less than zpt ensure proper 
vacuum-like behavior within the metal surface, where 
all electrons are explicit and thus no implicit screening 
should appear. 



B. Asymptotic behavior of electrostatic potential 

Unlike the standard Poisson equation, which has no 
unique solution for periodic systems because the zero of 
potential is an arbitrary constant, the modified Poisson- 
Boltzmann equation (jlOl) has a unique solution in pe- 
riodic systems. To establish this, we integrate the dif- 
ferential equation ([TO)) over the unit cell. The first term, 
which is the integral of an exact derivative, vanishes. The 
remaining terms then give the condition. 



K\r)q^{r)dV - — {Qr. 

£6 



(11) 



where Qn and Q n are the total number of electronic and 
nuclear charges in the cell, respectively. Any two (j){r) 
which differ by a constant C both can be valid solutions 
only if C / K^(r) dV = 0, so that we must have C = as 
long as K(r) is non-zero at any location in the unit cell. 
Thus, any amount of screening at any location in space in 
the calculation eliminates the usual indeterminacy of (j}{r) 
by an additive constant, thereby establishing an absolute 
reference for the zero of the potential. 

To establish the nature of this reference potential, we 
first note that deep in the fluid, far from the electronic 
system, the electron density approaches n(r) — and 
the dielectric constant and screening lengths attain their 
constant bulk values e(r) — >■ Cb and K{r) — ?> Kb. Under 
these conditions, the Green's function impulse response 
of PH]) to a unit point charge is 



exp (-Kfcr) 



(12) 



a Coulomb potential screened by the dielectric response 
of the solvent and exponentially screened by the presence 
of ions. Next, we rearrange ((TU)) so that the left-hand side 
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has the same impulse response as the bulk of the fluid but 
with a modified source term, 

-4^(PsolW+PextW) (13) 
where we have defined 

PextW^-^((fb-e(r))V2</,(r) 

- (Ve(r)) • (V0(r)) + e, (n^r) - nl) 0(r)). 

(14) 

The key step now is to note that all source terms clearly 
vanish in the bulk of the fluid where p^Qiir) — > 0, e(r) 
Cb = constant, and k{t) — >■ Kb- From the exponential 
decay of the Green's function (|12p and the vanishing of 
Pgpj+Pgxt bulk region of the fluid, we immediately 

conclude that 0(r) — deep in the fluid region, thereby 
establishing that the absolute reference of zero potential 
corresponds to the energy of an electron solvated deep in 
the fluid region. 



(100) surfaces of each of these metals were computed 
within a supercell representation with a distance of 10 
times the lattice constant of each metal (in all cases 
around 30 A) between surface slabs of thickness of 5 
atomic layers. For these initial calculations, we were 
very conservative in employing such large regions be- 
tween slabs to absolutely eliminate electrostatic super- 
cell image effects between slabs. We strongly suspect 
that smaller supercells can be used in the future. All cal- 
culations presented employ optimized^S norm-conserving 
Kleinman-Bylander pseudopotentials'^^ with single non- 
local projectors in the s, p, and d channels, a plane-wave 
cutoff energy of 30 H, and employ a 8 x 8 x 1 /c-point 
Monkhorst-Pack^ grids to sample the Brillouin zone. 

The JDFTx-calculated lattice constants of the bulk 
metals within both exchange-correlation approximations 
when using 8x8x8 fc-point grids are shown in Table H] 
Clearly, the LDA and GGA lattice constants both agree 
well with the experiment. Except where comparisons are 
specifically made with LDA results, all calculations in 
this work employ GGA for exchange and correlation. 

TABLE I: Cubic lattice constant (A) in conventional 
face-centered cubic unit cell 



C. Future Improvements 

While offering a computationally efficient and simple 
way to study electrochemistry, the approximate func- 
tional (O is highly simplified and possesses several lim- 
itations which the more rigorous approach of Section II 
overcomes by coupling of an explicit solvent model for 
Q,iq[Na{r)^^ to the electronic system through an ap- 
proach similar to the molecular pseudopotentials pro- 
posed by Kim et al^ Such limitations include the fact 
that because we employ a linearized Poisson-Boltzmann 
equation, we do not include the nonlinear dielectric re- 
sponse of the fiuid (which other approaches in the lit- 
erature to date also ignore^^i^i) or nonlinear saturation 
effects in the ionic concentrations, both of which become 
important for potentials greater than a few hundred mV. 
Despite these limitations, we remain encouraged by the 
promising results we obtain below for this simple func- 
tional and optimistic about the improvements that work- 
ing within a more rigorous framework would provide. 

V. ELECTRONIC STRUCTURE 
METHODOLOGY 

All calculations undertaken in this work and pre- 
sented in Section VI were all performed within the 
DFT-I--I- framework^^ as implemented in the open- 
source code JDFTxi^ They employed the local-density 
or generalized-gradient^^, approximations using a plane- 
wave basis within periodic boundary conditions. The 
specific materials under study in this paper were plat- 
inum, silver, copper, and gold. The (111), (110), and 



Metal LDA GGA Experimenl^i 



Pt 


3.93 


3.94 


3.92 


Cu 


3.55 


3.67 


3.61 


Ag 


4.07 


4.13 


4.09 


Au 


4.05 


4.14 


4.08 



VI. RESULTS 

To evaluate the promise of our approach, we begin by 
studying the fundamental behaviors of transition metal 
surfaces in equilibrium with an electrolyte environment 
as a function of applied potential. We find that even our 
initial highly simplified form of joint density-functional 
theory reproduces with surprising accuracy a wide range 
of fundamental physical phenomena related to electro- 
chemistry. Such transition metal systems, especially 
platinum, are of electrochemical interest as potential cat- 
alysts for both the oxygen reduction reaction (ORR) and 
th hydrogen evolution reaction (HER). Molecular dy- 
namics studies of the platinum system in solution, both 
at the classical-^ and ab initic^-^'^ levels, to date have 
not fully accounted for ionic screening in the electrolyte, 
which is essential to capturing the complex structure of 
the electrochemical double layer and the establishment 
of a consistent reference potential. 

For the initial exploratory studies presented in this 
manuscript, we focus on pristine surfaces without ad- 
sorbates in order to establish clearly the relationship be- 
tween theoretical and experimental quantities and to lay 
groundwork for future systematic comparison of potential 
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catalyst materials. Unless otherwise specified, we carry 
out our calculations with screening lengths of 3 A, cor- 
responding to monovalent ionic concentrations of 1.0 M. 
We employ these high concentrations because most elec- 
trochemical cells include a supporting electrolyte with 
high ionic concentration chosen to provide strong screen- 
ing while avoiding (to the extent possible) interaction 
with and adsorption on the electrode. Note that, because 
our present model includes only ionic concentrations and 
no other species-specific details about the ions in the 
electrolyte, our results correspond to neutral pH. Future 
work will readily explore pH and adsorption effects by in- 
cluding protons and other explicit ions in the electronic- 
structure portions of the calculation. One great advan- 
tage of the present theoretical approach is the ability to 
separate the role of the non- adsorbing ions in the sup- 
porting electrolyte from the role of the adsorbing ions 
that interact directly with the surface. 



A. Treatment of charged surfaces in periodic 
boundary conditions 

The application of voltage essential to the ab initio 
study of electrochemistry requires precise treatment of 
charged surfaces not accessible to common electronic 
structure approaches due to singularities associated with 
the Coulomb interaction. In the case of a vacuum envi- 
ronment, the electrostatic potential 4>{r) of even a neutral 
electrode approaches a physically indeterminate constant 
which varies with the choice of supercell. As is well- 
known, this difficulty compounds radically when a net 
charge is placed on the surface, resulting in a formally 
infinite average electrostatic potential in a periodically 
repeated system. By default, most electronic structure 
packages designed for use with periodic systems treat 
this singularity by setting the G = Q Fourier component 
of 0(r) to zero, equivalent to incorporating a uniform, 
neutralizing charge background throughout the region of 
the computation. This solution to the Coulomb infin- 
ity is not realistic in electrochemical applications where 
the actual compensating charge appears in the fiuid and 
should not be present in the interior of the electrode. 

Another option which has been employed in the elec- 
trochemical contexts^ is to include an oppositely charged 
counter-electrode located away from the working elec- 
trode in the vacuum region of the calculation. However, 
including an explicit density-functional electrode is of- 
ten computationally prohibitive as it requires doubling 
the number of electrons and atoms and requires a large 
supercell to prevent image interaction. Implicit inclu- 
sion of a counter-electrode through either coulomb trun- 
cation or an external charge distribution^^, requires an 
arbitrary choice of the distribution of external charges 
representing the counter-electrode, and such arbitrary 
choices may result in unphysical electrostatic potentials, 
even in the presense of a few explicit layers of neutral 
liquid molecules. One realistic choice is to employ De- 



bye screening as in Eq. [TO] This approach ensures that 
the long-range decay of 0(r) into the fiuid corresponds 
to the behavior of the actual physical system, that the 
fiuid response contains precisely the correct amount of 
compensating charge, and that the potential approaches 
an absolute reference, even in a periodic system. 

Another more explicit, and hence computationally ex- 
pensive, option employed in the electrochemical litera- 
ture is to add a few layers of explicit water molecules to 
the surface and then include explicit counter-ions (pro- 
tons) located in the first water layer^i^ This approach 
models some of the most important effects of the actual 
physical distribution of counter-ions, which really should 
contain both localized and diffuse components, by con- 
sidering only the first layer of localized ions. 

Figure Efa) contrasts the potential profiles resulting 
from the aforementioned approaches in actual calcula- 
tions of a Pt(lll) electrode surface. Figure[3Ka) displays 
the microscopic local electron potential energy function 
(0(z)) for a surface at applied voltage £ — —1.09^ vs 
PZC which corresponds to a charge of cr = —18 /iC/cm^. 
The screened electron potentials generated by solution 
of Eq. [TU] at two different ionic strengths (c = 1.0 M 
and c = 0.1 M) are compared to potential profiles for a 
similarly charged surface in vacuum, with the net charge 
in the system neutralized either by imposing a uniform 
background charge or by placing an oppositely charged 
counter electrode at one Debye screening length from the 
metal surface. The two charge-compensated vacuum cal- 
culations clearly do not correspond to the electrochemi- 
cal behavior, with far wider potential variations than ex- 
pected. Figure [3Kb) shows a detailed view of the macro- 
scopic electrostatic potential {^(z)) for the same charged 
surface (obtained by subtracting the microscopic electron 
potential of the neutral surface and switching the sign 
to reflect electrochemical convention) for the JDFT cal- 
culated charged surface at two different ionic strengths. 
The charge-compensated vacuum calculations would be 
off the scale of this figure, while the macroscopic elec- 
trostatic potential for the JDFT calculations obtains the 
value of the applied potential within the electrode and 
then approaches a well-established reference value of zero 
with the correct asymptotic behavior in the fluid region. 



B. Electrochemical double layer structure 

The Gouy-Chapman-Stern model, described in Section 
III(C), offers a well-known prediction of the structure 
of the electrochemical double layer, to which the poten- 
tials from our model correspond precisely. The electro- 
static potential profiles in the standard electrochemical 
picture include an initial, capacitor-like linear drop in 
(^(z)) due to the outer Helmholtz layer (the Stern re- 
gion), followed by a characteristic exponential decay to 
zero deep in the fluid (the diffuse Gouy-Chapman region). 
Our model naturally captures this behavior as a result of 
(a) the localization of the dielectric response and screen- 
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FIG. 3: Microscopic electron potential energies {4>{z)) 
and macroscopic electrostatic potentials {^{z)) 
averaged in planes for the Pt (111) surface as a function 
of distance z — zpt from the end of the metal surface: 
(a) {4>{z)) for surface with applied voltage f = — 1.09 V 
vs. PZC in vacuum (green dashed) and in monovalent 
electrolytes of c = 1.0 M (red) and c = 0.1 M (blue) 
where the dotted lines represent calculations with an 
explicit counter-electrode and the solid lines are JDFT 

calculations (b) close-up view of {^{z)) for JDFT 
calculations with c = 1.0 M (red) and c — 0.1 M (blue) 

and applied voltage £ = —1.09 V vs. PZC (almost 
indistinguishable in the previous plot) (c) Variation of 
{^{z)) in JDFT monovalent electrolyte of c = 1.0 M 
with £ = {-1.09, -0.55, 0.0, 0.55, 1.09} V vs. PZC. 



ing to the liquid region as described by Niq{r) through 
Eq. dS]) and (b) the separation between the fluid and re- 
gions of high explicit electron density n(r) through the 
definition of Niq{r) = Niq (n(r)) via Eq. ®. Both the 
Stern and Gouy-Chapman regions are clearly evident in 
Figures [3Jb,c). We find the dielectric constant transi- 
tion region appearing in Figure (Hb) , approximately the 
width of a water molecule, to be essential to the accurate 
reproduction of the double layer structure. The poten- 
tials for charged surfaces in Figures |31[b,c) first show a 
linear decay in the region A< z — zpt < A, corre- 
sponding to the "gap" between the end of the surface 
electron distribution (zpt) and the beginning of the fluid 
region, precisely the behavior we should expect in the 
Stern region. For a Pt(lll) surface at applied voltage 
-1.09 V vs. PZC, A = 0.6 A, but the width of this gap 
is voltage-dependent (as shown in Figure |4Kb)) and also 
varies with metal and crystal face. After the gap region, 
for A < z — zpt < A + J (where 7 = 0.6 as in^, the 
dielectric constant in Figure [H[b) changes rapidly from 
about 10 to the bulk value eb ~ 80, defining a transition 
region which ensures that no significant diffuse decay in 
the potential occurs until beyond the outer Helmholtz 
layer, thereby allowing proper formation of the diffuse 
Gouy-Chapman region for z — zpt > A -|- 7. We empha- 
size that we have not added these phenomena into our 
calculations a posteriori, but that they occur naturally 
as a consequence of our microscopic, albeit approximate, 
ab initio approach. 



C. Charging of surfaces with electrode potential 

To explore the effects of electrode potential on the sur- 
face charge and electronic structure. Figure Hl^a) shows 
the surface charge cr as a function of potential £ for a 
series of transition metal surfaces for an electrolyte of 
monovalent ionic strength c = 1.0 M, without adsorp- 
tion of ions to the surface. We find the average double 
layer capacitance of the Pt(lll) surface - the slope of the 
corresponding a — £ curve in Figure [DJ a) - to be C=19 
/iF/cm^, in excellent agreement with the experimental 
value of 20 fiF/cui'^M Indeed, we find that a significant 
fraction of our total capacitance is due to dielectric and 
screening effects in the fluid; this agreement again sup- 
ports our model for the electrolyte. The remainder is 
associated with the "quantum capacitance" or density of 
states CQQg of the surface slab in our supercell calcula- 
tions. 

Closer inspection of the charge versus potential data 
reveals that the slope is not quite constant as a func- 
tion of voltage. Indeed, taking the numerical derivatives 
of the curves in Figure HUa) yields values for the differ- 
ential capacitance that exhibit an approximately linear 
dependence on voltage. This voltage-dependence con- 
strasts with studies performed using a different technique 
to produce voltage-independent constant values for the 
capacitance*^, which not only was limited to producing 
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a constant value for the capacitance, but also required 
computationally demanding thermodynamic sampling to 
model the fluid. 

To understand the origin of the above voltage- 
dependence of the capacitance, we employ the series 
model for differential capacitance in ([3]) and Q , in which 
the total capacitance per unit area C is modeled as a 
series combination of the capacitance associated with 
the density of states of the metal, a Stern capacitance 
(Ca) across a gap of width A, and the (constant) Gouy- 
Chapnian capacitance associated with the inverse screen- 
ing length K. We can then extract the "gap" capacitance 
as 

— ^C^^=C ^^Cr)\)Q^——- (15) 

To verify that the voltage-dependence of this contribu- 
tion indeed correlates to changes in the gap associated 
with the Stern layer, we make an independent definition 
of the width of the gap as A = — ZMetai, where Zc 
represents the location where the presence of our model 
fluid becomes significant and ZMetai represents the loca- 
tion of the surface of the metal. Specifically, we define 
Zc as the point where the planar average of the inverse 
dielectric constant has fallen by half from its value in 
the electrode (as in Figure Sl^b)) since the polarization of 
the fluid becomes significant when (e~^(zc)) < 0.5. We 
determine z Metal from the covalent radii of the metal sur- 
face atoms, but note that the specific choice of ZMetai is 
unimportant in the analysis to follow. 

Figure IHc) correlates the inverse gap capacitance C^^ 
from the right hand side of Eq. [TS] with the values of A 
defined as in the previous paragraph. There is a strik- 
ing linear trend with a slope within about ten percent 
of , validating that the primary contribution to the 
voltage-dependence of the differential capacitance within 
this model is from changes in the gap between the fluid 
surface and where the dielectric screening begins. The 
ultimate origin of this effect within the present approxi- 
mation (in which the dielectric constant is determined by 
the electron density through Eq. [S]) can be traced to the 
increase in surface electrons with decreasing applied po- 
tential, which moves the location of the fluid transition 
further away from the metal surface. In fact, the ex- 
perimentally determined capacitance of Ft (111) due to 
only the double layer~ (after subtracting the effects of 
counter-ion adsorption) has a voltage-dependence quite 
similar to our prediction. Since the distance of closest 
approach of the fluid to the metal surface is determined 
by van der Waals interactions and the addition of more 
electrons could indeed strengthen the repulsion, the qual- 
itative voltage-dependence of the "double-layer" capaci- 
tance even at this simple level of approximation may in- 
deed be capturing some aspects of the underlying physics. 

The double-layer capacitance notwithstanding, in 
physical systems the total capacitance is dominated by 
the effects of adsorption of counter-ions, and so the quali- 
tative voltage-dependence of the capacitance at this sim- 



ple level of approximation has limited practical relevance. 
Nonetheless, it is an important feature of the electro- 
chemical interface for those modified Poisson-Boltzmann 
approaches in which the cavities are determined by con- 
tours of the electron density. Future work in this area 
could capture the "ion-adsorption" portion of the ca- 
pacitance either by including explicit counter-ions within 
the electronic structure portion of the calculation or by 
choosing a classical fluid functional that includes a mi- 
croscopic description of the counter-ions. 



D. Potentials of Zero Charge and reference to 
Standard Hydrogen Electrode 

To connect our potential scale (relative to an electron 
solvated in our model fluid) to a standard potential scale 
employed in the literature and to confirm the reliabil- 
ity of our model. Figures [5ja,b) show our ab initio pre- 
dictions for potentials of zero charge versus experimen- 
tal values relative to the standard hydrogen electrode 
(SHE) .4'' Within both the local density (LDA)i^ and gen- 
eralized gradient (GGA)^^ approximations to the elec- 
tronic exchange-correlation energy, we have calculated 
the potentials of zero charge for various crystalline sur- 
faces of Ag, Au, and Cu, three commonly studied metals. 
We performed a least-squares linear fit to the intercept 
of our data, leaving the slope fixed at unity. (Note that 
the experimental data for Cu in NaF electrolyte were 
not included in the fit, due to concerns discussed be- 
low.) The excellent agreement between our results (with 
a constant offset) and the experimental data indicates 
that joint density-functional theory accurately predicts 
trends in potentials of zero charge, and encourages us 
that it can establish oxidation and reduction potentials 
in the future. The improved agreement of GGA (rms er- 
ror: 0.058 V) over LDA (rms error: 0.108 V) underscores 
the importance of gradient-corrections to this type of sur- 
face calculation. 

The strong linear correlation with unit slope between 
the theoretical and experimental data in Figures [SJa,b) 
indicates that the simplified Poisson-Boltzmann ap- 
proach reproduces potentials of zero charge well rela- 
tive to some absolute reference. The single parameter 
in the fit for each of the two panels (namely, the ver- 
tical intercepts of each fit line) establishes the absolute 
relationship between our zero of potential (implicit in 
each set of theoretical results) and the zero of potential 
on the standard hydrogen-electrode scale (implicit in the 
experimental data). Specifically, we find that our zero 
of potential sits at -4.91 V relative to the SHE for LDA 
and -4.52 V relative to the SHE for GGA. Intriguingly, 
these values are close to the experimentally determined 
location of vacuum relative to the the standard hydrogen 
electrode reference (-4.44 V)^; in fact, the GGA refer- 
ence value is within a tenth of a volt. This apparent 
alignment is not altogether surprising due to the follow- 
ing argument: (1) our method measures the difference 
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FIG. 4: (a) Surface charge ct as a function of applied 
voltage £ for a series of transition metal surfaces in an 
electrolyte of monovalent ionic strength c = 1.0 M (b) 
Inverse dielectric constant e^^ as a function of distance 

from a Pt(lll) surface for multiple values of applied 
voltage (c) Inverse gap capacitance C^^ as a function of 
the distance from the metal surface at which the fluid 
begins A. The solid line indicates the best fit to the 
data with slope constrained to eg ^ 



in energy between an electron in the electrode and an 
electron solvated deep in our model electrolyte, so that 
our potentials of zero charge are measured relative to a 
solvated electron reference; (2) the energy of a solvated 
electron relative to vacuum within the presentlt consid- 
ered linearized Poisson-Boltzmann model is zero because 
this approximation includes only electrostatic effects; and 
(3) because the calculated potentials of zero charge in 
the figures are thus relative to vacuum, the difference be- 
tween our calculated results and the experimental results 
should represent the constant difference between the vac- 
uum and SHE references. 

Consideration of the breakdown of the potential of zero 
charge into physically meaningful quantities explains the 
difference between the LDA and GGA results and elu- 
cidates the apparent success of the rather simple mod- 
ified Poisson-Boltzmann approach in predicting PZC's. 
Transferring an electron from a metal surface to a refer- 
ence electrode requires, first, removal of the electron from 
the surface and, then, transport of the electron through 
the relevant interfacial layers of the liquid. The energy as- 
sociated with the fomer process is the work function, and 
the energy associated with the latter relates to the intrin- 
sic dipole of the liquid-metal interface. As is well-known, 
there is an approximately constant shift between the pre- 
dictions of the LDA and GGA exchange and correlation 
functionals for work functions of metals. In fact, Fall et 
al. report that GGA metal work functions are approxi- 
mately 0.4 V lower than the LDA work functions,-- corre- 
sponding well to the differences we find between the ver- 
tical intercepts of Figures [5l[a,b). Next, to aid considera- 
tion of the intrinsic dipole of the interface, Figure[SJc) ex- 
plicitly compares our predictions for work function with 
our predictions for potential of zero charge, including also 
the corresponding experimental data for both quantities. 
(To place all values on a consistent scale of potential, 
which we choose to be vacuum, we have added the ex- 
perimentally determined 4.44 V difference between SHE 
and vacuum to the experimental PZC's.) The data in 
Figure [5fc) suggest that the vacuum work functions are 
harder to predict than potentials of zero charge, possibly 
due to difficulty determining the value of the reference 
potential in the vacuum region, an issue not present in 
our fiuid calculations due to the screening in Eq. [TUl The 
figure also indicates an approximately constant shift from 
vacuum work function to potential of zero charge, sug- 
gestive of a roughly constant interfacial dipole for each 
of the metal surfaces. However, the shift is not exactly 
constant: both the experimental and theoretical data 
exhibit significant fluctuations (on the order of 0.1 V) 
in the shift between work function and PZC from one 
metal surface to another. Because the PZCs are deter- 
mined to within a significantly smaller level of fluctuation 
(0.06 V), these data indicate that the Poisson-Boltzmann 
model captures not merely a constant interfacial dipole, 
but also a significant fraction of the fluctuation in this 
dipole from surface to surface. 

We note that Tripkovic et al. have also calculated the 
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potentials of zero charge for transition metal surfacesJ^ 
However, that approach requires calculation of several 
layers of explicit water within the electronic structure 
portion of the calculation, and those authors find the 
resulting potentials of zero charge to be dependent on 
the exact structure chosen for the water layers. How- 
ever, while differing orientations of water molecules at 
the interface may result in significant local fluctuations 
in the instantaneous PZC, the experimentally measured 
potential of zero charge is a temporal and spatial thermo- 
dynamic average over all liquid electrolyte configurations 
rather than the value from any single configuration. Di- 
rect comparison to experimental potentials of zero charge 
therefore should involve calculation of a thermodynamic 
average. 

As a matter of principle, derivatives of the free en- 
ergy (which the JDFT framework provides directly) 
yield thermodynamic averages. Therefore, an exact free- 
energy functional would predict the exact, thermody- 
namically averaged potential of zero charge, and clas- 
sical liquid functionals, which capture more microscopic 
details of the equilibrium liquid configuration^^ than the 
present model, would be an ideal choice for future in- 
depth studies. Indeed, such functionals are capable of 
capturing the relevant electrostatic effects even when a 
single configuration of water molecules dominates the 
thermodynamic average. (In such cases, minimization 
of the free-energy functional results in localized site den- 
sities Na(r) representing the dominant liquid configura- 
tion.) Of course, in cases of actual charge-transfer reac- 
tions between the surface and the liquid, the (relatively 
few molecules) involved in the actual transfer must be 
included within the explicit electronic density-functional 
theory, whereas the other electrolyte molecules may still 
be handled accurately within the more computationally 
efficient liquid density-functional theory. 

There is also reason to be sanguine regarding the abil- 
ity of the modified Poisson-Boltzmann approximation 
pursued in this work to capture intcrfacial dipole effects. 
The macroscopic dielectric constant contained within the 
present model describes primarily the orientational polar- 
izability of water, so that the liquid bound charges result- 
ing from the minimization of the free energy should re- 
fiect the most dominant configurations of water molecules 
in the thermodynamic average, even if only a single con- 
figuration dominates. On the electrode side, the image 
charges corresponding to the bound charge also natu- 
rally appear, as a consequence of both the electrostatic 
coupling in our model and the metallic nature of the sur- 
face described within electronic density-functional the- 
ory. From an optimistic perspective, it is quite pos- 
sible that a significant portion of the electrostatics of 
the surface dipole would be captured even at the sim- 
plified level of a Poisson-Boltzmann description. Ulti- 
mately, quantification of how much of the effect is cap- 
tured may only be verified by comparison to experiment. 
For the systems so far considered, the excellent a pri- 
ori agreement between experimental measurements and 



our theoretical predictions indicates that the relevant ef- 
fects are indeed captured quite well. It appears that even 
a simple continuum model (which only accounts for the 
effects of bound charge at the interface and the corre- 
sponding image charges within the metal) can predict 
accurately key electrochemical observables such as po- 
tential of zero charge. Certainly, for more detailed fu- 
ture studies, we would recommend exploring the perfor- 
mance of more explicit functionals. However, the appar- 
ent accuracy and computational simplicity of the cur- 
rent Poisson-Boltzmann approach render it well-suited 
for high-throughput studies of electrochemical behavior 
as a function of electrode potential. 

As a further example of the utility of the Poisson- 
Boltzmann approach, the potential of zero charge calcu- 
lation for copper illustrates how this theory can be used 
as a highly controlled in-situ probe of electrochemical 
systems, with the ability to isolate physical effects which 
are not possible to separate in the experiment. Specifi- 
cally, in Figures [5] (a) and (b), for copper there are ex- 
perimental values for two different electrolytes, NaF and 
KCIO4, which are both claimed to be noninteracting with 
the metal surface4i Clearly, our theoretical values, which 
correspond to potentials of zero free charge without ad- 
sorption of or chemical reaction with ions from the elec- 
trolyte, agree more favorably with experimental data for 
the Cu surface in KCIO4 than for the NaF electrolyte. 
Our results suggest that future experimental exploration 
is warranted to investigate potential interactions between 
the NaF electrolyte and the copper surfaces, or into other 
possible causes of the discrepancy in potentials of zero 
charge. Perhaps some polycrystalline impurities caused 
the experimental potentials of zero charge of the sup- 
posed single-crystalline faces to become much more sim- 
ilar than our calculations and the KCIO4 data indicate 
they should be. Ab initio calculations offer an avenue to 
study each of these potential causes independently and 
to elucidate the mechanisms underlying the apparent ex- 
perimental disagreement. 

Finally, although potentials of zero charge are quite 
readily observed in experiments for less reactive met- 
als such as silver and gold, measurement of the poten- 
tial of zero charge for platinum can be difficult because 
platinum is easily contaminated by adsorbates. For this 
reason, more convoluted methods are employed to de- 
termine an experimental value for the potential of zero 
charge for platinum. For instance, one may turn to ultra- 
high vacuum methods, where, by definition, no molecules 
are adsorbed on the surface, and one may then attempt 
to estimate the effect of the solution on the potential 
of zero free charge4i Alternately, one may employ cyclic 
voltammograms to estimate the charge due to adsorbates 
and then extrapolate the potential of zero free charge.— 
Our ab initio method, however, gives the values for non- 
contaminated potentials of zero free charge directly, pro- 
vided we establish the relation of our zero reference po- 
tential relative to that of the standard hydrogen elec- 
trode, which we have done above. 
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FIG. 5: Comparisons of ab initio predictions and 
experimental datai^ for potentials of zero free charge 
(PZCs) and vacuum work functions: (a) ab initio LDA 
predictions versus experimental PZC relative to SHE, 
(b) ab initio GGA predictions versus experimental PZC 

relative to SHE, (c) ab initio GGA vacuum work 
functions (solid line with squares) and PZCs (solid line 
with circles), experimental work functions (dotted line) 
and PZCs (dashed line) versus vacuum for the same 
series of surfaces. Best linear fits with unit slope (dark 
diagonal solid lines in (a) and (b)). 



For uncontaminated platinum, our method yields the 
potentials of zero free charge shown in Table [Hi Com- 
pared to other references in the literature, the best agree- 
ment with our results is from an experiment which ex- 
trapolates the potential of zero charge from ultra-high 
vacuum, eliminating the effects of unknown adsorbates 
on the clean surface.™ As with the results for Cu, the sig- 
nificantly better agreement of our calculations with this 
latter experimental approach suggests that perhaps fu- 
ture experiments which measure potential of zero charge 
should reconsider the effect of possible contaminants 
when extrapolating values for the potential of zero free 
charge. 

TABLE II: Platinum Potentials of Zero Free Charge 
(V vs SHE) 



(110) (100) (111) 
LDA 0.31 0.70 0.71 
GGA 0.40 0.79 0.82 



VII. CONCLUSION 

In this work, we extend joint density-functional the- 
ory (JDFT) - which combines liquid and electronic free- 
energy functionals into a single variational principle for 
a solvated quantum system - to include ionic liquids. 
We describe the theoretical innovations and technical de- 
tails required to implement this framework for study of 
the volt age- dependence of surface systems within stan- 
dard electronic structure software. We establish a con- 
nection to the fundamental electrochemistry of metallic 
surfaces, accurately predicting not only potentials of zero 
charge for a number of crystalline surfaces for various 
metals but also an independent value for the standard 
hydrogen electrode relative to vacuum. Furthermore, we 
show how future innovations in free energy functionals 
could lead to even more accurate predictions, demon- 
strating the promise of the joint density-functional ap- 
proach to predict experimental observables and capture 
subtle electrochemical behavior without the computa- 
tional complexity required by molecular dynamics simu- 
lations. These advantages render joint density-functional 
theory an ideal choice for high-throughput screening cal- 
culations and other applications in materials design. 

We have built extensively upon the framework of 
joint density-functional theory in the implicit solvent 
approximation^ extending it to include charged ions in 
a liquid electrolyte. Beginning with an implicit model for 
the fluid density Niq (r) in terms of the electronic density 
of the surface, Niq{r) — Niq{n{r)), we include an ionic 
screening length tied to the fluid density Niq{r) in the 
same way as done in previously successful models for the 
dielectric constant. We also solve a previously unrecog- 
nized difficulty by including model core electron densities 
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within the surface to prevent artificial penetration of hq- 
uid density into the ionic cores, which lack electrons in 
typical pseudopotential treatments of the solid. Inclusion 
of this ionic screening allows us to provide a consistent 
zero reference of potential and to resolve many difficulties 
associated with net charges in periodic supercell calcula- 
tions, thereby enabling study of electrochemical behavior 
as a function of applied voltage. 

With the framework to include electrode potential 
within joint density-functional theory calculations thus 
in place, we then establish clear connections between 
microscopic computables and experimental observables. 
We identify the electronic chemical potential of density- 
functional theory calculations with the applied voltage 
in electrochemical cells, and thereby extract a numerical 
value of 4.52 V (within the GGA exchange-correlation 
functional) for the value of the standard hydrogen elec- 
trode relative to vacuum, which compares quite favor- 
ably to the best-accepted experimental value of 4.44 V.'*'* 
We also show that joint density-functional theory repro- 
duces, a priori, the subtle voltage-dependent behaviors 
expected for a microscopic electrostatic potential within 
the Gouy-Chapman-Stern model and we extract poten- 
tials of zero free charge for a series of metals commonly 
studied in electrochemical contexts, often finding agree- 
ment with experimental values to within hundredths of 
volts. 

This qualitatively correct prediction of electrochemi- 
cal behavior and encouraging agreement with experiment 
demonstrate the capabilities of even a simple approxi- 
mation within the joint density-functional theory frame- 
work, and we expect future improvements to the free- 
energy functional to be able to describe more complex 
electrochemical phenomena. Future work should also 
generalize the approximate functional to include nonlin- 
ear saturation effects in ionic screening within the cur- 
rent modified Poisson-Boltzmann approach, with an ap- 
proach along the lines of other works. In electrochem- 
ical experiments, the differential capacitance of charged 
metal surfaces often exhibits a minimum at the poten- 
tial of zero charge^S' (not seen in the linear continuum 
theory) , and more advanced theories including such non- 
linear effects should be able to capture this more subtle 
behavior. Additionally, recent developments in classical 
density functionals for liquid water— now can be imple- 
mented to study electrochemical systems. Such classi- 
cal density functionals can be extended to include re- 
alistic descriptions of ions and are capable of capturing 
other essential behaviors of electrolyte fluids, including 
features in the ion-ion and ion-water correlation func- 
tions due to differences in the structure of the anion and 
the cation. Finally, in systems where electrochemical 
charge-transfer reactions are important or where chemi- 
cal bonds of the fluid molecules are expected to break, the 
relatively few reactant molecules should be treated within 
the explicit electronic structure portion of the calcula- 
tion, with the remaining vast majority of non-reacting 
molecules handled within the more computationally effi- 



cient liquid density-functional theory. 

With advances such as those described above, joint 
density-functional theory holds promise to become a use- 
ful and versatile complement to the toolbox of currently 
available techniques for first principles study of electro- 
chemistry. Unlike ah initio molecular dynamics (or any 
other theory involving explicit water molecules), this 
computationally efficient theory is not prohibitive for 
larger system sizes. In fact, as the system size grows, 
the fraction of calculation time spent solving the modified 
Poisson-Boltzmann equation actually decreases, meaning 
that for larger systems, the calculation is only nominally 
more expensive than calculations of the corresponding 
systems carried out in a vacuum environment. Also, 
because thermodynamic integration is not required, the 
joint density-functional theory approach yields equilib- 
rium properties directly and has a clear advantage over 
molecular dynamics simulations for calculation of free 
energies. Immediate applications include the study of 
molecules on metallic electrode surfaces as a function of 
applied potential and prediction of the basic properties 
of novel catalyst and catalyst support materials. These 
calculations could inform future materials design by of- 
fering an opportunity to screen novel complex oxides and 
intermetallic materials in the presence of the true elec- 
trochemical environment, thereby elucidating the funda- 
mental physical processes underlying fuel cells and liquid- 
phase Graetzel solar cells. 
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Appendix A: Implementation within standard 
electronic structure software 



Here we consider the issues which arise when im- 
plementing the above framework within a pre-existing 
electronic-structure code. We will focus on software oper- 
ating within the pseudopotential framework as this tech- 
nique is commonly used for surface calculations. 

Such pseudopotential calculations, for computational 
efficiency, include the nuclei and core electrons together 
as a unit and describe their combined effects on the 
valence electron system through effective, "pseudo"- 
potentials. Two subtleties now arise. First, because the 
pseudopotentials describe the long-range electrostatic in- 
teraction between the ionic cores and the electrons, the 
screening of the long-range Coulomb part of the pseu- 
dopotentials by the electrolyte environment through (fTO)) 
must be handled properly. Second, because the calcu- 
lated (valence) electron density n(r) in the atomic core 
regions tends to be relatively low in pseudopotential cal- 
culations, our definition of the liquid density Niq{n) as a 
local function of the local electron density n{r) through 
([9|) can lead to the unphysical presence of liquid within 
the atomic cores if precautions are not taken. 
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As a matter of notation specific to this appendix, we 
separate conceptually the valence electron density n^{r), 
calculated directly with the Kohn-Sham orbitals, from 
the missing contribution nc{r, {Ri}) due to the core elec- 
trons, which clearly varies explicitly with the locations 
of the centers of the ions. By including both electron 
and ionic (now, actually, valence-electron and ionic core) 
source terms, the new energy functional ([7]) naturally 
provides electrolyte screening of all of the relevant fields. 
The functional ([7]) then becomes 

A[n{r),4>{r)] = AtxcKW] 

+ ScK(r),n,(r,{i?,}),</)(r)] 

+ Ups[nv{r),{Zj,Rj}], (Al) 

where Atxc['>T'v{'>')] is the Kohn-Sham single-particle ki- 
netic plus exchange-correlation energy, and 

ScKW,ne(r,{i?,}),0(r)] = 

dM0W {nv{r)~N{r,{Zi,Ri})) 

ebK^{ny{r) +nc{r,{Ri})) 2 



Stt 



iHr)r} (A2) 



represents the previously described electrostatic 
contributions to the total energy functional (with 
N{r,{Zi,Ri}) = J2iZiS^^Hr - Ri) representing point 
charges with the ionic valences Zj), and the term 

Ups[nv{r),Zi,Ri] = 

/ (r^AV^^Pir-Ri)^ n^ir) d\, (A3) 



with 

AV^Pir - Rj) = Vp^Pir - Ri) + ZiG{r - Ri), (A4) 

represents the non-Coulombic components of the pseu- 
dopotential VpP [r), with G{r) = l/\r\ being the Coulom- 
bic Green's function associated with a unit point charge 
in free space. Note that here and henceforth in this 
work Zj refers to the charges of the ionic pseudopoten- 
tial cores and not the nuclear atomic numbers. Finally, 
Appendix iBl details how, for practical numerical reasons, 
we work not with mathematical point charges but rather 
with narrow charge distributions which we can resolve 
numerically. 

The derivative of the functional ^[71(7-)] with respect to 
ny{r) at a point r (which gives the local effective Kohn- 
Sham potential used for the electronic wavefunction min- 
imization) must also be adjusted to include the new di- 



electric response and ionic screening terms, 
d d 



dny{r 



+ Y,AVp^P{r-Rj) 
I 



(A5) 



Finally, Hellman-Feynman calculation of the forces on 
the atoms requires care because of the dependence of 
both the ionic core density N{r,{Zi,Ri}) and the model 
core electron density nc{r,{Ri}) on the ionic positions 
Rj. The final result is 

Vfl,,A = Vr, Ec + Vi^, Ups, (A6) 

where 

Vfl, Ups = J {Vr, AVpP\r - i?,)) n,(r) d^r 

X Vfl,nc(r, {i?/}) 

d^r(l){r)ZiVRj'^^\r~Ri) 



(A7) 



While some of these derivatives have simpler analyti- 
cal forms than those described above, these forms render 
much simpler the numerical representation of the rel- 
evant quantities, particularly the Dirac-delta functions 
and pseudopotentials, as described in Appendix IB] 



Appendix B: Numerical Details 

The system-dependent modified Poisson-Boltzmann 
equation which appears in our calculations does not have 
a direct analytic solution in either Fourier or real space 
and, thus, requires a numerical solution such that we can- 
not employ analytic Dirac S functions to represent the 
ion-core charges. Instead, in our numerical calculations, 
we employ "smoothed" ion-core charge densities. 



iV(^)(^ {Zi,Ri}) = Y,ZiS^''\r~Ri) 



(Bl) 



where (S^"^) is a normalized, isotropic three-dimensional 
Gaussian of width a containing a single unit change. To 
the extent that the smoothed distributions do not overlap 
with each other or the fluid regions, the replacement of 
the point charges with these distributions will not affect 
the "Ewald" energy among the charged atomic cores or 
the dielectric screening effects. In practice, we find good 
numerical solutions by employing a relatively narrow 
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width determined by the spatial resolution of the calcula- 
tion, so that a corresponds to a distance of 1.40 points on 
the Fourier grid. (Specifically, this work employs a plane- 
wave energy cutoff of 30 H, so that a = 0.168 A.) This 
choice of parameters ensures that there is little overlap 
among the Gaussians and between the Gaussians and the 
fluid, so that that this replacement has negligible effect 
on the screened interaction among the atomic cores. 

These smoothed distributions, however, do overlap 
with the valence electrons, an effect which we must com- 
pensate. We compensate this local effect exactly by re- 
placing the point-charge response G(r) in the modifica- 
tion of the pseudopotential (|A4[) with the response cor- 
responding to the smoothed densities, 

g(.)(,),erf(M/V2^)^ 
\r\ 

We also represent the core-electron densities in Eq. 
that prevent fluid penetration into the atomic cores, and 
whose form is thus not critical, with Gaussian distribu- 
tions 

n,(r,{i?,}) = C^5('-c)(r-i?,). (B3) 
/ 

Alternately, for this density, one could use the core- 
electron density from the partial core correction from an 
appropriately designed pseudopotential. Because the role 
of ric in our framework is simply to prevent penetration 
of fluid into the ionic cores, the precise values of the norm 
C and width rc are not critical. We find that the choice 
C = 0.3 A~^, rc = 0.2 Aworks well for this purpose for 
all the species in our calculations. 

Finally, with the above definitions in place, we have 
taken care to make all replacements S^{r) — >■ ^("'^(r), 
G{r) -> G'^'^^r) and nc(r,{i?/}) = C Y.iS^''^Hr - Ri) 
in the appropriate places in the expressions for the total 
free energy, in the functional derivatives appearing in the 



effective Kohn-Sham potential (jA5l) . and in the expres- 
sions for the Hellman-Feynman forces on the atoms (|A6p . 
These substitutions complete the numerical specification 
of the functionals employed in our calculations. 

We find that standard electronic structure methods 
work well with our functionals. The one equation whose 
solution requires new algorithms is the modified Poisson- 
Boltzmann equation, which, unlike the standard Pois- 
son equation, does not have a direct analytic solution in 
Fourier space. To solve this equation, we have, however, 
found a simple to implement, yet highly efficient precon- 
ditioned conjugate gradient algorithm. 

The portions of the functional which depend on the 
potential field (j){r) appear in Ec in Eq. (jA2[) . which is 
a quadratic functional whose maximum corresponds to 
the solution of the modified Poisson-Boltzmann equation, 
Eq. (1101) . and whose quadratic kernel is 

Q = (V-eV-efcK2)/(47r). (B4) 

We chose to solve this equation in Fourier space, where 
the diagonal elements of this kernel have a very simple 
approximate form, which leads to the diagonal precondi- 
tioner, 

K{G) = {eG^ + ebK'^)-\ (B5) 

where e = ^ / ^{^)<i^i" a-nd = ^ J K^{r)(Pr are the 
average values of these parameters over the unit cell. 
This diagonal preconditioner completely ignores the spa- 
tial variation of the dielectric constant. A more effec- 
tive preconditioner which may be obtained by building 
in this inhomogeneity - is calculated by first multiply- 
ing by y/jK{G)) (the square root of the diagonal pre- 
conditioner) in Fourier space, transforming to real space 
and dividing by e(r), then returning to Fourier space and 
again multiplying by y/XK{G)). This inhomogeneous 
preconditioner requires more time to evaluate for a single 
iteration than the diagonal preconditioner, but reduces 
the total number of iterations required significantly. 



